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Abstract 

Neurotrophins are a diverse class of structurally related proteins, essential for neuronal development, survival, 
plasticity and regeneration. They are characterized by major family members, such as the nerve growth factors 
(NGF), brain-derived neurotrophic factors (BDNF) and neurotrophin-3 (NT-3), which have been demonstrated here to 
lack coding sequence variations and follow the regime of negative selection, highlighting their extremely important 
conserved role in vertebrate homeostasis. However, in stark contrast, venom NGF secreted as part of the chemical 
arsenal of the venomous advanced snake family Elapidae (and to a lesser extent Viperidae) have characteristics 
consistent with the typical accelerated molecular evolution of venom components. This includes a rapid rate of 
diversification under the significant influence of positive-selection, with the majority of positively-selected sites found 
in the secreted p-polypeptide chain (74%) and on the molecular surface of the protein (92%), while the core structural 
and functional residues remain highly constrained. Such focal mutagenesis generates active residues on the toxin 
molecular surface, which are capable of interacting with novel biological targets in prey to induce a myriad of 
pharmacological effects. We propose that caenophidian NGFs could participate in prey-envenoming by causing a 
massive release of chemical mediators from mast cells to mount inflammatory reactions and increase vascular 
permeability, thereby aiding the spread of other toxins and/or by acting as proapoptotic factors. Despite their 
presence in reptilian venom having been known for over 60 years, this is the first evidence that venom-secreted NGF 
follows the molecular evolutionary pattern of other venom components, and thus likely participates in prey- 
envenomation. 
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Introduction 

Venom, a complex biochemical cocktail of biologically active 
components, such as proteins, peptides, amino acids, 
neurotransmitters and polyamines, has underpinned the 
diversification and evolutionary success of several animal 
lineages [1]. This key evolutionary innovation is employed by a 
plethora of animals for predation, competitor deterrence and 
defence [2-5]. The scientific consensus is that venom 
components originate via toxin recruitment events, as part of 
which physiological protein-encoding genes are duplicated and 
the new copies are selectively expressed in the venom gland 
[5-15]. Over the years, our understanding of the origin and 
diversification of snake venoms has greatly increased, largely 
due to advances in transcriptomics and proteomics [16-22]. 
However, the precise role of certain proteins, which are 
secreted as part of the biochemical venom arsenal, still 
remains to be elucidated. Nerve growth factor (NGF), a key 
member of the neurotrophin family, is one such class of protein 
whose presence in snake venoms has been intriguing. Since 
its discovery in the late 1950s , NGF has been reported from 
the venoms of various caenophidian (advanced) snakes, 
including members of the front-fanged elapid and viperid 
families as well as from venomous lizards [23-29], but its 
function and relative importance in snake venoms remains 
unknown [8]. 

Neurotrophins represent a family of structurally related 
proteins, crucial for neuronal development, survival, death, 
regeneration and plasticity. According to the classical 
neurotrophic hypothesis, neurotrophins are produced in limiting 
amounts and the survival of the innervating neurons is 
dependent on winning the competition for sufficient quantities 
of these factors [30,31]. Neurotrophins contain gene family 
members such as nerve growth factor (NGF), brain-derived 
neurotrophic factor (BDNF), neurotrophin-3 (NT-3) and 
neurotrophin-4/5 (NT-4/5) [32-35], all of which function by 
interacting with the p75 neurotrophin receptor (NTR) in their 
proneurotrophin forms and various structurally related 
tropomyosin-related kinase or Tyrosine kinase receptors (Trk) 
in their active cleaved form [36,37]. While NGF specifically 
activates TrkA [38], BDNF and NT-4/5 interact with TrkB 
[36,39]. NT-3 primarily interacts with TrkC and is unique in also 
being capable of weakly binding to both TrkA and TrkB [36]. 

Neurotrophins have been extensively studied not only 
because they are perceived as one of the primary factors 
responsible for the complexity of vertebrate nervous systems, 
but also because of their involvement in cognition and memory. 
Knockout of the genes encoding NGF, BDNF and NT-3 genes 
is fatal in mice, highlighting the importance of these proteins for 
survival and normal neuronal development. Not surprisingly, 
abnormalities associated with the production of neurotrophins 
have been linked with neuropathies and neurodegenerative 
disorders. 

In order to investigate the role of NGF in the venom of 
Toxicofera reptiles [24], we have investigated the molecular 
evolution of these proteins in reptilian (turtles; squamates: 
Laterata, Scinciformata, Gekkota; Toxicofera lizards: 
Anguimorpha and Iguania; Henophidia snakes; advanced 



snakes: Elapidae, Viperidae and 'non-front-fanged' advanced 
snakes) and mammalian lineages, by employing sophisticated 
protein and codon-level selection assessments. We further 
compare the molecular evolution of NGF with the other major 
members of the neurotrophin family, namely BDNF and NT-3, 
in a wide array of reptilian and mammalian lineages. Molecular 
evolution analyses conducted on a dataset of 1183 nucleotide 
sequences revealed that these genes have remained largely 
unchanged since their origin over 300 million years (the split 
between mammals and reptiles: www.timetree.org) due to the 
extremely important functions they play in vertebrate 
homeostasis. 

Results 

Bayesian and maximum-likelihood analyses of NGF, BDNF 
and NT-3 genes retrieved trees with the same topology 
(Figures 1-3; Figure S1-S4), which were in concordance with 
the earlier reported phylogenies of neurotrophins [40]. 

One ratio model (ORM), the simplest of the lineage-specific 
models, computed w of less than 0.50, therefore indicating an 
influence of negative selection on NGF, BDNF and NT-3 genes 
in all of the mammalian and reptilian lineages examined 
(Tables S1-S3). This highly conservative model can only detect 
positive selection when the ui ratio averaged over all the sites 
along the lineages in a phylogenetic tree is significantly greater 
than one. Nevertheless, the computed w of 0.86 for the 
Elapidae NGF highlights a relatively greater accumulation of 
variation in this lineage (Table S1). 

To detect episodic diversifying selection, which only affects 
certain sites in the protein, we employed site-specific models 
(Tables 1-3; Tables S1-S3). Like ORM, site model 8 (M8) also 
indicated a strong influence of negative selection on various 
reptilian and mammalian NGF, BDNF and NT-3 genes (Tables 
1-3; Tables S1-S3; Figures 4 and 5). However, Elapidae NGF 
was found to be rapidly evolving (to = 0.95). The Bayes 
Empirical Bayes (BEB) approach in M8 identified as many as 
19 positively-selected sites (13% of the total sites with w of 
4.09) in Elapidae NGF, indicating the strong influence of 
positive selection on this lineage (Table 1 and Table S1). 
Although, M8 failed to detect variation in the Viperidae (w = 
0.77) and 'non-front-fanged' advanced snake (w = 0.58) NGF 
lineages, as many as 9% (w = 2.40) and 8% (w = 2.50) of the 
total codon sites were detected as rapidly diversifying in the 
respective lineages. Single Likelihood Ancestral Counting 
(SLAC), Fixed Effects Likelihood (FEL), Random Effects 
Likelihood (REL), Mixed Effects Model of Evolution (MEME), 
Fast, Unconstrained Bayesian AppRoximation (FUBAR) and 
integrative approach, conclusively supported these findings 
and highlighted the complete lack of variations in NGF, BDNF 
and NT-3 genes in various reptilian and mammalian lineages 
examined (Tables 1-3). However, elapid venom NGF was 
found to be rapidly evolving under the influence of positive 
selection (integrative analyses: 34 positively selected sites). 
Viperidae NGF was also found to accumulate relatively greater 
variations (integrative analyses: 8 positively selected sites) in 
comparison to the non-venomous reptilian and mammalian 
NGF lineages. Analysis of NGF from the venom-glands of 'non- 
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Figure 1. Bayesian molecular phylogeny of nerve growth factors (NGF). Branches with the Bayesian posterior probability 
(B.P.P) of less than 0.85 are highlighted in grey (remaining in colours). Site model 8 (M8) computed oj values for respective 
lineages are presented. The number of positively selected sites (PP > 0.95) detected by M8's Bayes-Empirical Bayes (BEB) 
approach in Elapidae lineage is also indicated. Elapid sequences representing putative duplicate genes are indicated with red labels 
[NFF: "non-front-fanged" advanced snakes; SCI: Scinciformata]. 

doi: 10.1371/joumal.pone.0081827.g001 



front-fanged' advanced snakes and venomous lizards 
(Anguimorpha and Iguania) was hampered by the scarcity of 
sequences. Despite this, the integrative approach was able to 
identify five positively selected sites in the 'non-front-fanged' 
advanced snake lineages. 



We further evaluated selection pressures along the elapid 
and viperid NGF lineages using two-ratio model (which 
assesses selection pressures only across lineages) and the 
branch-site test A (which assesses selection pressures across 
the sites and along the lineages). The two-ratio model failed 
(p>0.05) to detect positive selection in elapid (w=0.85) and 
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Figure 2. Bayesian molecular phylogeny of brain-derived neurotrophic factors (BDNF). Branches with the Bayesian posterior 
probability (B.P.P) of less than 0.85 are highlighted in grey (remaining in colours). Site model 8 (M8) computed oj values for 
respective lineages are presented. 

doi: 10.1371/joumal.pone.0081827.g002 



viperid NGF (w=0.82). However, branch-site test A indicated a 
strong influence of positive selection on elapid NGF [w = 3.90; 
4.5% positively selected (PS) sites; significant at 0.001 after 
Bonferroni correction] and to a lesser extent on Viperidae NGF 
(w = 2.01 and 8.1% of PS-sites; significant at 0.001 after 
Bonferroni correction: Table S4). We further employed the 
clade model c approach to simultaneously compute and 
compare the w values of various Toxicofera NGF lineages 
(Table 1). Clade model analyses indicated the influence of 



positive selection in shaping the evolution of Elapidae NGF 
(w=2.38), while the remaining lineages were found to be 
constrained by negative selection (but again with the caveat 
that fewer 'non-front-fanged' advanced snake and Toxicofera 
lizard venom-gland specific NGF sequences were available; 
Table 1). 

To derive further support for the sites detected as positively 
selected by the nucleotide analyses of elapid NGF, we 
employed a complementary amino acid-level approach 
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Figure 3. Bayesian molecular phylogeny of neurotrophln-3 (NT-3). Branches with the Bayesian posterior probability (B.P.P) of 
less than 0.85 are highlighted in grey (remaining in colours). Site model 8 (M8) computed w values for respective lineages are 
presented. 

doi: 10.1371/journal.pone.0081827.g003 



implemented in TreeSAAP (Table 4). Since a mutated amino 
acid can share similar/identical biochemical and structural 
properties with the ancestral amino acid being replaced, not all 
non-synonymous mutations may affect the structure and 
function of a protein. Hence, the evaluation of the influence of 
selection on the biochemical and structural properties of 
proteins provides further insights about the strength of a 



mutation in affecting the fitness of an organism. With this 
combined approach, we were able to identify 13 positively 
selected sites (68% of positively selected sites) in elapid NGF 
(Table 4). Mapping of these mutations onto the three- 
dimensional structure of NGF revealed that a majority of 
detected hypermutable sites (74%) were concentrated in the 
secreted region (B-polypeptide chain) of Elapidae NGF (Figure 
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Table 1. Molecular evolution of Nerve Growth Factor (NGF). 









SLAC a 


FEL b 


REL C 


FUBAR d 


Integrative 


MEME e 


M8 


M2a 


Clade 


Elapidae 




I0>1* 


2 


13 


28 


22 


34 




19 


15 






tu<ia 


12 


20 


11 


18 


22 


21 


(14+5) 


(13+2) 


2.38 




10= 


0.90 




1.93 








0.95 


0.94 




Viperidae 




I0>1* 


0 


0 


2 


0 


8 




0 


0 






I0<18 


3 


13 


0 


7 


13 


6 






1.03 




10= 


0.65 




0.88 








0.77 


0.77 




NFF 




I0>1 f 


0 


1 


3 


2 


5 




0 


0 






I0<1S 


1 


2 


0 


3 


3 


2 






0.57 




co= 


0.59 




0.67 








0.58 


0.58 




Henophidia 




I0>1 f 


0 


0 


0 


0 


1 




0 


0 






I0<1S 


3 


12 


1 


13 


13 


1 






0.0001 




co= 


0.30 




0.43 








0.30 


0.30 




Iguania 




I0>1* 


0 


0 


0 


0 


4 




0 


0 






io<ia 


35 


55 


73 


83 


87 


4 






0.20 




co= 


0.26 




0.76 








0.25 


0.29 




Anguimorpha 




I0>1* 


0 


0 


0 


0 


1 




0 


0 






L0<1S 


4 


13 


0 


15 


15 


1 






0.20 




10= 


0.33 




0.48 








0.33 


0.36 




Gekkota 




UJ>1 f 


0 


0 


3 


1 


5 




0 


0 






I0<1S 


1 


15 


0 


16 


16 


1 










co= 


0.24 




0.34 








0.25 


0.25 




Scinciformata 




OJ>1 f 


0 


0 


0 


0 


1 




0 


0 






I0<1S 


8 


18 


All 


18 


19 


1 










co= 


0.22 




0.26 








0.22 


0.25 




Laterata 




I0>1* 


0 


0 


0 


0 


0 




0 


0 






to<1fl 


5 


12 


0 


14 


15 


0 










10= 


0.26 




0.36 








0.26 


0.27 




Turtles 




I0>1* 


0 


0 


0 


1 


1 




0 


0 






to<1S 


9 


13 


0 


14 


14 


0 










co= 


0.33 




0.38 








0.32 


0.34 




Mammals 




I0>1 f 


0 


1 


0 


0 


4 




0 


0 






I0<1S 


24 


54 


All 


86 


86 


3 










10= 


0.21 




0.22 








0.17 


0.11 





S5). Mutation mapping also revealed that a greater proportion identified a significantly large proportion of sites in elapid NGF 

of amino acid residues (92%; remaining sites couldn't be as evolving under the influence of positive selection, while the 

assigned to exposed/buried class) were exposed to the remaining reptilian and mammalian NGF, BDNF and NT-3 

surrounding medium, indicating that these Elapid venom lineages were depicted as constrained by extreme negative 

proteins evolve through focal mutagenesis (Table 4). selection pressures (Figure S6-S8). It is noteworthy that a 

Evolutionary fingerprint analyses, which is based on the small fraction of sites in all snake lineages (Caenophidia and 

probability distribution of site-to-site synonymous (a or dS) and Henophidia) seem to have accumulated a greater proportion of 

nonsynonymous ((3 or dN) substitution rates in an alignment, non-synonymous to synonymous mutations, although these 
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Table 1 (continued). 



a: Single Likelihood Ancestor Counting 

b: Fixed-effects likelihood 

c: Random-effects likelihood 

d: Fast, Unconstrained Bayesian AppRoximation 

Integrative: Sites detected in common by SLAC, FEL, REL, FUBAR and MEME 

e: Sites detected as experiencing episodic diversifying selection (0.05 significance) by the Mixed Effects Model Evolution (MEME)M8: Positively-selected sites detected 
using the Bayes Empirical Bayes approach implemented in M8. Sites detected at 0.99 and 0.95 significance are indicated in the parenthesisM2a: Positively-selected sites 
detected using the Bayes Empirical Bayes approach implemented in M2a. Sites detected at 0.99 and 0.95 significance are indicated in the parenthesisClade: Omega 
computed by the clade model 

f: Number of positively selected sites at 0.05 significance (for SLAC, FEL) or 50 Bayes factor (for REL) / number of sites under pervasive diversifying selection at the 
posterior probability >0.9 (FUBAR) 

g: Number of negatively selected sites at 0.05 significance (for SLAC, FEL) or 50 Bayes factor (for REL) / number of sites under pervasive purifying selection at the posterior 
probability >0.9 (FUBAR) 

u>: mean dN/dSNFF: "non-front-fanged" advanced snakes 
doi: 10.1371/joumal.pone.0081827.t001 



observations lacked significant statistical support (indicated by 
the density and compactness of the circles). Branch-site REL 
identified four, three and two branches (p<0.05), respectively in 
elapid, viperid and 'non-front-fanged' advanced snake NGFs as 
evolving under the influence of episodic diversifying selection. 
This highlights the episodic nature of natural selection on the 
evolution of advanced snake venom-specific nerve growth 
factors (Figure S9). This test failed to identify positive 
diversifying selection in other non-venomous reptilian or 
mammalian NGF and any BDNF and NT-3 gene lineages, 
including those of the advanced snakes. 

Discussion 

Distinct rates and trajectories of neurotrophin evolution 

Evidence provided by various analyses (codeml site, branch, 
branch-site and clade-specific models: M8, M2a, M3, M0, two- 
ratio model, branch-site test A, clade model c; HyPhy: SLAC, 
FEL, REL, MEME, FUBAR, integrative analyses, branch-site 
REL; TreeSAAP: amino acid-level selection assessment and 
the evolutionary fingerprint analysis) not only revealed the 
strong influence of positive diversifying selection on Elapidae 
venom NGF (and to a lesser extent on the venom NGF of other 
advanced snake lineages, particularly the Viperidae NGF), but 
also highlighted the role of purifying selection pressures in 
shaping the molecular evolution of various non-venomous 
reptilian and mammalian NGF lineages (Tables 1-3; Figures 
1-5; Tables S1-S4; Figures 4 and 5; Figures S6-S9). The 
accumulation of variation in advanced snake NGFs, in a 
fashion similar to other venom proteins [5,7,8] is suggestive of 
a role in envenomation and prey-capture. As envenomation is 
their primary method of prey subjugation, elapid and viperid 
snakes rely on the swift and potent action of venom for 
foraging. Hence, most elapid and viperid venom components 
experience extreme selection pressures, likely as a result of a 
co-evolutionary arms race with their prey [5,7]. Although NGF 
constitutes a very small proportion (1-5 mg/g of venom) of the 
total venom in most snakes [29,41^3], the venom of snakes of 



the genus Oxyuranus (Taipans) has been shown to contain 
larger quantities of NGF [44^6]. NGF seems to make up 
nearly 0.5% (w/w) of the total venom injected by these snakes 
[44-46]. Thus, NGF could be secreted in elevated amounts by 
several species of caenophidian snakes, which could be 
enough to cause toxicity in the prey. 'Non-front-fanged' 
advanced snakes have been shown to possess several rapidly 
evolving venom components that participate in prey 
envenoming [47-52]. Selection analyses in the present study 
identify several residues and branches in both the 'non-front- 
fanged' and viperid advanced snake NGF lineages as evolving 
under positive Darwinian selection (Table 1: integrative 
approach; Table S4: branch-site and clade models; Figure S6, 
S9b and S9c). However, due to the scarcity of NGF sequences 
from both these lineages, the extent to which NGF participates 
in their envenoming remains elusive. 

The aforementioned analyses also demonstrated the 
influence of negative selection on the evolution of all reptilian 
and mammalian BDNF and NT-3 gene lineages examined in 
this study. The lack of variation in these genes despite their 
origin over 300 million years ago is probably due to their 
significance in neuronal development, maintenance and 
survival. NGF, BDNF and NT-3 deficiencies have been 
correlated with severe peripheral neuropathy and death of the 
organism [53]. 

Focal mutagenesis of venom nerve growth factors 

The synthesis and secretion of venom proteins is an 
energetically expensive process [54-56]. Hence, mutations that 
disrupt the structure/function of proteins are filtered out of the 
population by negative selection pressures, favouring the 
conservation of catalytic and structurally important core 
residues. As a result, most mutations in proteins can be found 
concentrated in structurally and/or functionally unimportant 
regions. Moreover, the accumulation of point mutations under 
the influence of positive Darwinian selection in certain regions 
of the toxin, such as the molecular surface (a phenomenon we 
refer to as focal mutagenesis), may confer adaptive 
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Table 2. Molecular evolution of Brain-derived Neurotrophic Factors (BDNF). 









SLAC a 


FEL b 


REL C 


FUBAR d 


Integrative 


MEME e 


M8 


M2a 


Elapidae 




I0>1* 


0 


0 


0 


0 


0 




0 


0 




U)<19 


1 


3 


0 


3 


3 


0 








10= 


0.09 




0.25 








0.09 


0.09 


Viperidae 




I0>1* 


0 


0 


0 


0 


0 




0 


0 




io<ia 


1 


3 


13 


4 


13 


0 








10= 


0.10 












0.11 


0.11 


NFF 




I0>1 f 


0 


0 


0 


0 


0 




0 


0 




I0<1S 


5 


19 


14 


22 


22 


0 








co= 


0.08 




0.35 








0.06 


0.06 


Typhlopoidea 




I0>1 f 


0 


0 


0 


0 


0 




0 


0 




I0<1S 


27 


47 


0 


93 


93 


0 








10= 


0.07 




0.11 








0.07 


0.07 


Boidae 




I0>1* 


0 


0 


1 


1 


1 




0 


0 




to<1S 


3 


10 


3 


10 


10 


0 








co= 


0.0 




0.28 








0.09 


0.09 


Pythonidae 




I0>1* 


0 


0 


0 


1 


1 




0 


0 




L0<1S 


0 


2 


11 


7 


11 


0 








10= 


0.21 




0.23 








0.15 


0.15 


Iguania 




I0>1 f 


0 


0 


0 


0 


0 




0 


0 




U)<19 


24 


56 


111 


99 


24 


1 








10= 


0.05 








0.05 




0.05 


0.05 


Anguimorpha 




I0>1 f 


0 


0 


2 


0 


0 




0 


0 




io<ia 


9 


22 


9 


23 


9 


0 








10= 


0.08 




0.15 




0.08 




0.06 


0.06 


Crocodiles 




I0>1* 


0 


0 


0 


0 


0 




0 


0 




L0<1S 


0 


8 


20 


1 


0 


0 








10= 


0.09 








0.09 




0.12 


0.12 


Turtles 




I0>1* 


0 


0 


4 


0 


5 




0 


0 




io<ia 


6 


18 


2 


18 


18 


1 








C0= 


0.10 




0.14 








0.05 


0.06 


Mammals 




I0>1 f 


0 


0 


2 


0 


4 




0 


0 




io<ia 


48 


72 


51 


113 


72 


2 








10= 


0.05 












0.05 


0.06 



significance. We and others have postulated that predatory area (ASA) ratio for each residue indicated that 92% of the 

venom components undergo focal mutagenesis that results in positively selected residues in these proteins are located on the 

the generation of a diversity of novel amino acids (and their molecular surface (Table 4; remaining sites couldn't be 

side-chains) on the molecular surface, which could non- conclusively assigned to exposed/buried class). Moreover, a 

specifically interact with prey cells and cause a myriad of large proportion of these positively selected sites (74%) were 

pharmacological effects [11,52,57-61]. Mapping of found concentrated in the p-chain of Elapidae NGF, which is 

hypermutable sites onto the three-dimensional homology the only secreted region of the mature toxin and hence the only 

model of elapid NGFs and the calculation of accessible surface region that is likely to experience a coevolutionary arms race 
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Table 2 (continued). 



a: Single Likelihood Ancestor Counting 

b: Fixed-effects likelihood 

c: Random-effects likelihood 

d: Fast, Unconstrained Bayesian AppRoximation 

Integrative: Sites detected in common by SLAC, FEL, REL, FUBAR and MEME 

e: Sites detected as experiencing episodic diversifying selection (0.05 significance) by the Mixed Effects Model Evolution (MEME) 

M8: Positively-selected sites detected using the Bayes Empirical Bayes approach implemented in M8. Sites detected at 0.99 and 0.95 significance are indicated in the 
parenthesis 

M2a: Positively-selected sites detected using the Bayes Empirical Bayes approach implemented in M2a. Sites detected at 0.99 and 0.95 significance are indicated in the 
parenthesis 

f: Number of positively selected sites at 0.05 significance (for SLAC, FEL) or 50 Bayes factor (for REL) / number of sites under pervasive diversifying selection at the 
posterior probability >0.9 (FUBAR) 

g: Number of negatively selected sites at 0.05 significance (for SLAC, FEL) or 50 Bayes factor (for REL) / number of sites under pervasive purifying selection at the posterior 
probability >0.9 (FUBAR) 

uj: mean dN/dSNFF: "non-front-fanged" advanced snakes 
doi: 10.1371/journal.pone.0081827.t002 



(Figure S5). Thus, it is highly likely that elapid NGFs follow the 
regime of focal mutagenesis and favour the accumulation of 
hypermutable sites in specific regions of the secreted toxin that 
are able to influence the fitness of the animal. Rapidly evolving 
caenophidian NGFs with a diversity of amino acids on their 
molecular surface could not only induce a plethora of 
pharmacological effects by non-specifically interacting with 
novel receptors of the prey, but could further aid in evading the 
prey immune response [62]. 

Putative role of nerve growth factors in venom 

Nerve growth factors are known to inhibit venom 
metalloproteinase-dependent proteolysis [63]. Hence, the 
presence of large amounts of NGF, relative to 
metalloproteinases in viper venom [64,65] has been viewed as 
an indication of their role in prevention of venom autolysis [63]. 
However, additional plausible mechanisms of action exist for 
both the direct and indirect participation of rapidly evolving 
caenophidian NGF in prey envenoming, potentially resulting in 
an increase in the potency of the toxicoferan venom arsenal. 
The injection of large amounts of NGF into the lymphatic 
system of prey would result in NGF overdose, triggering a 
variety of reactions. Nerve growth factors are also known to 
cause apoptosis of cells lacking the expression of TrkA 
receptors [66-68] and thus the caenophidian venom NGFs 
could act as cytotoxic proapoptotic factors, after being 
introduced into tissues where they are not typically found. 
Venom components like hyaluronidase have been shown to 
degrade hyaluronic acid in the extra-cellular matrix and function 
as spreading factors, facilitating the easier diffusion of other 
venom components [69]. NGF, which is known to trigger the 
release of granules containing histamines, serotonins and other 
chemical mediators from mastocytes [70,71], could perform a 
similar ancillary venom function by facilitating the efficient 
absorption of other venom components. The degranulation of 
mast cells that are located primarily in perivascular spaces, 
often close to neurons and blood vessels, has been correlated 



with increased vascular permeability and even neurogenic 
inflammations [72]. While controlled degranulation aids in 
orchestrating and mounting acute inflammatory reactions, 
massive releases can be associated with anaphylaxis, 
bronchoconstriction (suffocation) and vasodilation [73]. Over 
the years, both glycosylated and non-glycosylated forms of 
NGFs have been reported from snake venoms [44,74-76]. 
Glycosylation has been proposed to prolong the circulation 
time of serum proteins by increasing their stability [77]. 
Interestingly, Asn 23, which has been postulated as a putative 
glycosylation site in mature venom NGF [78,79], is highly 
conserved in the entire squamate lineage (except in a few 
elapid snakes: genus Naja, Notechis scutatus, Pseudechis 
australis, Pseudechis porphyriacus, etc.), while being 
completely absent from mammalian and turtle NGFs (Figure 
S5). These observations are not surprising since mammalian 
NGFs are not known to undergo glycosylation. Snake venom 
NGFs could undergo posttranslational glycosylation in order to 
circulate longer in the bloodstream of the prey and spread 
further throughout the prey animal's system. Thus, 
caenophidian NGFs could directly and/or indirectly participate 
in prey-envenoming through a number of plausible 
mechanisms. 

Putative duplication of nerve growth factors in Elapidae 

Genes encoding venom proteins are known to evolve via the 
birth-and-death model of evolution where new genes are 
created by repeated duplication events, and are subsequently 
either maintained in the genome, deleted or become non- 
functional pseudogenes [7,8]. The molecular phylogeny of 
NGFs clearly indicates the presence of two independent NGF 
genes in Naja sputatrix, with very strong support for the node 
separating them (915 bootstrap and 1.0 Bayesian posterior 
probability: Figure 1 ). The phylogenetic placement of these two 
gene homologues suggests that the gene duplication event 
responsible for their origin occurred in an early elapid ancestor 
(Figure 1) - we might therefore expect to find additional copies 
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Table 3. Molecular evolution of of Neurotrophin 3 (NT-3). 









SLAC a 


FEL b 


REL C 


FUBAR d 


Integrative 


MEME e 


M8 


M2a 


Elapidae 




u)>1* 


0 


0 


0 


0 


0 




0 


0 




w<1S 


1 


4 


0 


4 


4 


0 








co= 


0.19 




0.21 








0.18 


0.19 


Viperidae 




I0>1 f 


0 


0 


0 


0 


0 




0 


0 




U)<1S 


0 


5 


0 


3 


5 


0 








U)= 


0.30 




0.49 








0.28 


0.30 


NFF 




I0>1 f 


0 


2 


5 


1 


8 




0 


0 




uxia 


43 


61 


79 


82 


95 


2 










0.22 




0.35 








0.20 


0.27 


Typhlopoidea 




I0>1* 


1 


2 


8 


2 


10 




0 


0 




w<ia 


48 


69 


49 


106 


106 


2 








U)= 


0.24 




0.49 








0.23 


0.26 


Boidae 




I0>1* 


0 


0 


0 


1 


3 




0 


0 




0)<1S 


7 


23 


10 


26 


26 


2 








co= 


0.24 




0.51 








0.24 


0.28 


Scinciformata 




I0>1 f 


0 


0 


0 


0 


0 




0 


0 




uxia 


7 


27 


75 


30 


75 


0 










0.13 












0.12 


0.13 


Iguania 




U)>1 f 


0 


2 


4 


3 


7 




0 


0 




uxia 


27 


41 


23 


58 


58 


3 








0J= 


0.25 




0.43 








0.23 


0.26 


Anguimorpha 




U)>1 f 


0 


0 


0 


0 


1 




0 


0 




axis 


7 


20 


All 


23 


23 


1 










0.24 




0.32 








0.21 


0.24 


Crocodiles 




I0>1 f 


0 


0 


0 


0 


0 




0 


0 




uxia 


0 


1 


all 


1 


1 


0 










0.08 




0.08 








0.07 


0.07 


Turtles 




U)>1 f 


0 


0 


0 


0 


0 




0 


0 




uxia 


3 


4 


all 


5 


5 


0 










0.19 




0.21 








0.17 


0.17 


Mammals 




U)>1 f 


0 


0 


1 


0 


2 




0 


0 




uxia 


77 


124 


134 


167 


172 


1 










0.06 












0.05 


0.07 



of NGF in other elapid snakes in the future. Although one might be exclusively associated with a role in prey envenoming. We 

initially assume that the original copy of NGF would continue therefore propose that co-expression of an NGF gene occurs in 

the ancestral role in homeostasis, thereby releasing the novel different tissue types, similar to that recently proposed for other 

duplicate gene to evolve a function specific to venom, in N. toxin families [81]. Specifically, we postulate that despite the 

sputatrix both NGF genes were found to be expressed in the elevated mutation rate in the Elapidae NGF, the ancestral copy 

venom gland [80]. However, considering how important NGF of the gene has retained its biochemical functions by 

and other neurotrophins appear to be for the survival of an undergoing focal mutagenesis and continues to ensure 

organism [53], it is highly unlikely that both gene copies would homeostasis. A greater accumulation of mutations on the 
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Table 3 (continued). 



a: Single Likelihood Ancestor Counting 

b: Fixed -effects likelihood 

c: Random-effects likelihood 

d: Fast, Unconstrained Bayesian AppRoximation 

Integrative: Sites detected in common by SLAC, FEL, REL, FUBAR and MEME 

e: Sites detected as experiencing episodic diversifying selection (0.05 significance) by the Mixed Effects Model Evolution (MEME) 

M8: Positively-selected sites detected using the Bayes Empirical Bayes approach implemented in M8. Sites detected at 0.99 and 0.95 significance are indicated in the 
parenthesis 

M2a: Positively-selected sites detected using the Bayes Empirical Bayes approach implemented in M2a. Sites detected at 0.99 and 0.95 significance are indicated in the 
parenthesis 

f: Number of positively selected sites at 0.05 significance (for SLAC, FEL) or 50 Bayes factor (for REL) / number of sites under pervasive diversifying selection at the 
posterior probability >0.9 (FUBAR) 

g: Number of negatively selected sites at 0.05 significance (for SLAC, FEL) or 50 Bayes factor (for REL) / number of sites under pervasive purifying selection at the posterior 
probability >0.9 (FUBAR) 

uj: mean dN/dSNFF: "non-front-fanged" advanced snakes 
doi: 10.1371/journal.pone.0081827.t003 



Elapidae 

03=0.95 
P.S=19 


Viperidae 

03=0.77 
P.S=0 


'Non-front-fanged' 

05=0.58 
P.S=0 


Henophidia 

C0=0.30 
P.S=0 

Gekkota 


Anguimorpha 

(0=0.33 
P.S=0 

Scinciformata 


Iguania 

03=0.25 
P.S=0 

Laterata 



03=0.25 
P.S=0 



03=0.22 
P.S=0 



03=0.26 
P.S=0 



Mammals 




Turtles _ 

cted 03=0.32 
__ P.S=0 



| Positively selected 

mnm 

Variable Average Conserved 

Evolutionary Conservation 

Figure 4. Molecular evolution of nerve growth factors. Three-dimensional homology models of nerve growth factors, depicting 
the locations of positively selected sites (in red). The oj values (Model 8) for the respective taxa along with the number of positively 
selected sites (PP > 0.95, Bayes-Empirical Bayes approach) are indicated. 

doi: 10.1371/journal.pone.0081827.g004 
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Elapidae 



Non-front-fan gecP 
Advanced Snakes 




Pythonidae 



Boidae 



Typhlopidae 





Anguimorpha 




Iguanh 



Mammals 





(0-0.05 
P.S=0 




Turtle 



~ (0=0.32 
P.S=0 




fl | 2 | 3 I 4 | 5 | 6 | 7 |TB 

variable Average conserved 

Evolutionary Conservation 

Figure 5. Molecular evolution of brain-derived neurotrophic factors. Three-dimensional homology models of brain-derived 
neurotrophic factors (BDNF), depicting the molecular evolution of these non-venomous homologues. The a> values (Model 8) for the 
respective taxa (PP > 0.95, Bayes-Empirical Bayes approach) are indicated. 

doi: 10.1371/journal.pone.0081827.g005 



molecular surface and the loops of the protein ensures the 
preservation of structurally and functionally important core 
residues. Clearly, a large number of residues across the entire 
vertebrate NGF lineage appear to be highly conserved (Figure 
S5). Moreover, it has been demonstrated on numerous 
occasions that snake venom NGFs exhibit biological activities 
similar to their mammalian counterparts [44,74]. This would 
seemingly free the duplicate elapid NGF gene to participate 
exclusively in prey envenoming. In contrast to the Elapidae 
NGF, only one copy was found in the venom glands of the 
Anguimorpha lizard Abronia graminea, which was identical to 
the one recovered from genomic sequencing for use in 
taxonomy studies [27]. Moreover, despite extensive BLAST 
searches using sequence templates from a wide array of 
mammalian and reptilian NGFs against the genomes of Anolis 
carolinensis (Iguania) and Python molurus (Henophidia), we 
were only able to retrieve a single copy of the NGF gene. 
Hence, we further speculate that the ancestral NGF gene is co- 



expressed in physiological tissues and the venom-gland, while 
the new gene duplicates are more likely to be tissue- 
specifically expressed in the venom gland alone. Experimental 
investigation in the future regarding the copy number of NGF 
should reveal if the NGF duplication occurred in an early elapid 
ancestral lineage or at the base of the Toxicofera phylogenetic 
tree. 

Note on the usage of NGF in the phylogenetic 
reconstructions of Toxicofera reptiles 

NGF is a commonly used gene in reptilian phylogenetic 
studies [82-85]. However, the detection of gene duplication 
events and the accelerated evolution of the advanced snake 
NGFs indicate that these genes should not be employed in 
phylogenetic inferences of Toxicofera reptiles, as these factors 
can undermine such analyses [86]. 
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Table 4. Nucleotide and complementary protein-level selection assessment of Elapidae Nerve Growth Factors (NGF). 



Codon 


Amino Aacid 


M2a a 


M8 b 


Property 0 


Magnitude 01 


ASA 


57 


R 


4.241±0.515 


3.894±0.510 


- 


- 


- 






(0.994) 


(0.998)* 








72 


A 


4.260±0.453 


3.899±0.493 












(0.999)** 


(1.0)** 








116 


T 


4.227±0.553 


3.892+0.513 


Ra> a m 


6, 7 


■ 






(0.990)* 


(0.997)** 








123 


A 


4.243±0.509 


3.895+0.506 


Rat Cf/n 


6, 7 


■ 






(0.994) 


( 0.998)* 








125 




4.132±0.769 


3.857±0.601 


- 


HI 


IB 






(0.961)* 


(0.986)* 








127 


N 


4.262±0.448 


3.900+0.492 


Rai a m 


6, 7 


■ 






(1.0)** 


(1.0)** 








133 


Q 


4.262±0.446 


3.900+0.491 


Rat a m 


6, 7 


■ 






(1.0)** 


(1.0)** 








144 


D 


4.258±0.459 


3.899±0.495 


Rat a m 


6, 7 


69.5 






(0.999)** 


(1.0)** 






Exposed 


150 


T 


4.252±0.481 


3.897+0.500 


°m 


6 


77.5 






(0.997)** 


(0.999)** 






Exposed 


155 


R 


4.146±0.740 


3.867+0.576 




6 


87.0 






(0.965)* 


(0.990)* 






Exposed 


158 


V 


4.247±0.498 


3.896±0.504 


- 




76.9 






(0.995)** 


(0.999)** 






Exposed 


163 


E 


4.257±0.464 


3.898±0.496 




6 


65.9 






(0.998)* 


(0.999)** 






Exposed 


167 


L 


4.024±0.940 


3.824±0.670 




6 


65.0 






(0.929) 


(0.976)* 






Exposed 


170 


E 


4.260±0.453 


3.900±0.493 






73.0 






(0.999)** 


(1.0)** 






Exposed 


182 


R 


3.951 ±1.036 


3.789±0.736 






57.3 






(0.907) 


(0.965)* 






Exposed 


207 


Q 


4.255±0.471 


3.898+0.499 


Ra 


6 


32.2 






(0.998)** 


0.999** 






NA 


211 


R 


4.260±0.455 


3.899+0.493 


Ra 


6 


53.1 






(0.999)* 


(1.0)** 






Exposed 


219 


Q 


4.262±0.446 


3.900+0.491 


Ra 


6 


73.8 






(1.0)** 


(1.0)** 






Exposed 


240 


D 


4.004±0.970 


3.804+0.712 


Ra 


6 


66.9 






(0.923) 


(0.969)* 






Exposed 


167 


L 


4.024±0.940 


3.824±0.670 


°m 


6 


65.0 






(0.929) 


(0.976)* 






Exposed 


170 


E 


4.260±0.453 


3.900±0.493 






73.0 






(0.999)** 


(1.0)** 






Exposed 


182 


R 


3.951 ±1.036 


3.789±0.736 






57.3 






(0.907) 


(0.965)* 






Exposed 


207 


Q 


4.255±0.471 


3.898+0.499 


Ra 


6 


32.2 






(0.998)** 


0.999** 






NA 


211 


R 


4.260±0.455 


3.899+0.493 


Ra 


6 


53.1 






(0.999)* 


(1.0)** 






Exposed 


219 


Q 


4.262±0.446 


3.900+0.491 


Ra 


6 


73.8 






(1.0)** 


(1.0)** 






Exposed 


240 


D 


4.004±0.970 


3.804+0.712 


Ra 


6 


66.9 






(0.923) 


(0.969)* 






Exposed 
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Table 4 (continued). 



Amino-acid property symbols used: Power to be in the middle of a-helix (a m ), Solvent accessible reduction ratio (R a ) 
PAML 

a, b: Bayes Empirical Bayes (BEB) posterior probability and post-mean omega (indicated in brackets) for the sites detected as positively selected by the site models M2a 
and M8, respectively. Sites detected as positively selected at 0.95 and 0.99 posterior probability by the Bayes Empirical Bayes approach of M8 are represented by * and **, 
respectively. 
TreeSAAP 

c: amino acid property experiencing positive diversifying selection 
d: magnitude of selection on the amino acid property 
ASA: Accessible surface area. 

Note: Codon sites with significant support from both nucleotide and protein-level selection analyses are highlighted in bold, 
doi: 10.1371/journal.pone.0081827.t004 



Conclusions 

In this study, we show that BDNF and NT-3, two major 
members of the neurotrophin family, lack variations completely 
and evolve under the regime of negative selection, indicating 
their extreme importance in the maintenance of homeostasis. 
In complete contrast, this study points for the first time towards 
venom-specific NGF evolving in a manner typical of toxic 
proteins and consequently suggests a hitherto unappreciated 
role in envenomation. We highlight an exquisite mechanism of 
venom evolution where focal mutagenesis has facilitated the 
transformation of a key regulatory protein into a toxin, without 
hampering its participation in homeostasis. Structure-function 
studies of mutant NGF may reveal a suite of novel activities, 
some of which may be of use in drug design and development. 

Materials and Methods 

Phylogenetic Analyses 

Phylogenetic analyses were performed to allow the 
reconstruction of the molecular evolutionary history of 
vertebrate NGF, BDNF and NT-3 genes. Nucleotide sequences 
were downloaded from the National Center for Biotechnology 
Information database (NCBI: http://www.ncbi.nlm.nih.gov/) and 
a dataset of 1183 sequences (NGF: 308; BDNF: 434; NT-3: 
441) was compiled. Accession numbers of all sequences 
analysed in this study are available in Tables S5.1-S5.3. 
Unpublished sequences are available in the File S1. Resultant 
sequence sets were aligned using MUSCLE [87]. The 
phylogenetic relationships were determined using Bayesian 
and maximum-likelihood approaches. MrBayes version 3.2 [88] 
was used for Bayesian inference. Tree searches were run 
using four Markov chains for a minimum of 10 million 
generations, sampling every 100th tree. The log likelihood 
score of each saved tree was plotted against the number of 
generations to establish the point at which the log-likelihood 
scores of the analyses reached their asymptote. 25% of the 
total trees sampled were conservatively discarded as burnin. 
The posterior probabilities for clades were established by 
constructing a majority rule consensus tree for all trees 
generated after the completion of the burnin. The analyses 
were repeated three times to make sure that the trees 



generated were not clustered around local optima. An optimal 
maximum likelihood phylogenetic tree was obtained using 
PhyML 3.0 [89] and node support was evaluated with 1,000 
bootstrapping replicates. 

Test for Recombination 

To overcome the effects of recombination on the 
phylogenetic and evolutionary interpretations [90], we 
employed Single Breakpoint algorithm implemented in the 
HyPhy package and assessed recombination on all the toxin 
forms examined in this study [91,92]. When potential 
breakpoints were detected using the small sample Akaike 
information Criterion (AlCc), the sequences were 
compartmentalized before conducting selection analyses, so as 
to allow recombining units to have distinct phylogeny. 

Selection Analyses 

We evaluated the influence of natural selection on various 
members of vertebrate neurotrophins using maximum- 
likelihood models [93,94] implemented in CODEML of the 
PAML [95]. In order to detect evolutionary selection pressures 
acting upon individual lineages, we employed the two-ratio 
model as well as the optimized branch-site test [96,97]. The 
two-ratio model evaluates selection across the lineages alone, 
while the branch-site model allows oj to vary across the sites of 
the protein and along the branches in the tree. The latter is 
known for its reasonable power and accuracy to detect short 
bursts of episodic adaptations [97]. However, both the two-ratio 
and branch-site models require the foreground (lineages under 
positive selection) and background lineages (lineages lacking 
positive selection) to be defined a priori. Such predefined 
biological hypotheses are often unavailable and it becomes 
difficult to define the foreground lineages. Therefore, we 
treated each clade being compared as foreground branches 
alternatively and tested multiple hypotheses. A likelihood-ratio 
test was then conducted by comparing the model that allows oj 
to be greater than 1 in the foreground branch, with a null model 
that does not (oj constrained 1). It has been suggested that 
while implementing multiple hypotheses using branch and 
branch-site models, it is necessary to control the family-wise 
error rate (FWER or Type I error) [97]. We used Bonferroni 
correction to account for such errors. It uses a/n as the 
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significance level to test each hypothesis; where 'a' is the 
significance level and 'n' is the number of independent true null 
hypotheses. We further utilized branch-site REL implemented 
in HyPhy to identify lineages affected by positive selection 
pressures. Unlike the aforementioned lineage-specific models, 
branch-site REL does not require the identification of 
foreground and background branches a priori. 

Because lineage-specific models either assume a single oj 
for the entire tree or assess the effects of selection only along 
the branches, they often fail to identify regions in proteins that 
might be affected by episodic selection pressures and 
ultimately underestimate the strength of selection. Hence, we 
employed site-specific models which estimate positive 
selection statistically as a non-synonymous-to-synonymous 
nucleotide-substitution rate ratio (oj) significantly greater than 
1 . We compared likelihood values for three pairs of models with 
different assumed oj distributions as no a priori expectation 
exists for the same: MO (constant oj rates across all sites) 
versus M3 (allows the oj to vary across sites within 'ri discrete 
categories, n > 3); M1a (a model of neutral evolution) where all 
sites are assumed to be either under negative (uj < 1) or 
neutral selection (oj = 1) versus M2a (a model of positive 
selection) which in addition to the site classes mentioned for 
M1a, assumes a third category of sites; sites with oj >1 
(positive selection) and M7 (P) versus M8 (P and uj), and 
models that mirror the evolutionary constraints of M1 and M2 
but assume that oj values are drawn from a p distribution [98]. 
Only if the alternative models (M3, M2a and M8: allow sites 
with to >1) show a better fit in Likelihood Ratio Test (LRT) 
relative to their null models (MO, M1a and M7: do not allow 
sites oj >1), are their results considered significant. LRT is 
estimated as twice the difference in maximum likelihood values 
between nested models and compared with the x 2 distribution 
with the appropriate degree of freedom (i.e., the difference in 
the number of parameters between the two models). The 
Bayes empirical Bayes (BEB) approach [97] was used to 
identify amino acids under positive selection by calculating the 
posterior probabilities that a particular amino acid belongs to a 
given selection class (neutral, conserved or highly variable). 
Sites with greater posterior probability (PP > 95%) of belonging 
to the 'oj > 1 class' were inferred to be positively selected. 

SLAC, FEL, REL and FUBAR [99,100] implemented in 
HyPhy [101] were employed to provide additional support to 
the aforementioned analyses and to detect sites evolving under 
the influence of positive and negative selection. MEME [102] 
was also used to detect episodic diversifying selection. Further 
support for the results of the nucleotide-level selection 
analyses was obtained and the radicalness of mutations were 
assessed using a complementary protein-level approach 
implemented in TreeSAAP [103]. 

Direct comparison of uj values computed from different 
datasets can be misleading, as they can have different 
proportion of sites under selection. Hence, we assessed the 
selection pressures shaping the evolution of NGF along various 
toxicoferan reptilian lineages by employing clade model 
analyses implemented in codeml and simultaneously estimated 
oj values [104]. The significance of the analysis was tested by 
comparing the likelihood of this model with that of model M1a. 



To clearly depict the proportion of sites under different regimes 
of selection, an evolutionary fingerprint analysis was carried out 
using the evolutionary selection distance (ESD) algorithm 
implemented in datamonkey [91,105,106]. Evolutionary 
fingerprint analysis fits a versatile general discrete bivariate 
model of site-to-site variation in selection pressures and 
comprises a description of the number of selective classes, the 
dN/dS (oj) ratio for each class and the assignment of sites to 
classes [91,106]. 

Structural Analyses 

To depict the natural selection pressures influencing the 
evolution of various neurotrophins, we mapped the sites under 
positive selection on the homology models created using Phyre 
2 webserver [107]. Pymol 1.3 [108] was used to visualize and 
generate the images of homology models. The consurf 
webserver [109] was used for mapping the evolutionary 
selection pressures on the three-dimensional homology 
models. GETAREA [110] was used to calculate the ASA ratio 
or the solvent exposure of amino-acid side chains. It uses the 
atom co-ordinates of the PDB file and indicates if a residue is 
buried or exposed to the surrounding medium by comparing 
the ratio between side-chain ASA and the "random coil" values 
per residue. An amino-acid is considered to be buried if ASA is 
less than 20% and exposed if ASA is > 50%. 

Supporting Information 

Table S1. (S1.1 - S1.11) Details of selection analyses of 
nerve growth factors (NGF). a: dn/ds (weighted average), b: 
Significance of the model in comparison with the null model, c: 
Number of sites with oj > 1 under the Bayes empirical Bayes 
approach with a posterior probability (PP) more than or equal 
to 0.99 and 0.95. * Models which allow oj > 1 . 
(PDF) 

Table S2. (S2.1 - S2.11) Details of selection analyses of 
brain-derived neurotrophic factors (BDNF). a: dn/ds 
(weighted average), b: Significance of the model in comparison 
with the null model, c: Number of sites with oj > 1 under the 
Bayes empirical Bayes approach with a posterior probability 
(PP) more than or equal to 0.99 and 0.95. * Models which allow 
oj > 1. 
(PDF) 

Table S3. (S3.1-S3.11) Details of selection analyses of 
neurotrophin-3 (NT3). a: dn/ds (weighted average), b: 
Significance of the model in comparison with the null model, c: 
Number of sites with oj > 1 under the Bayes empirical Bayes 
approach with a posterior probability (PP) more than or equal 
to 0.99 and 0.95. * Models which allow oj > 1 . 
(PDF) 

Table S4. Lineage-specific selection analyses of nerve 
growth factors (NGF). a: dn/ds (weighted average), b: 
Significance of the model in comparison with the null model. * 
Significant after Bonferroni correction. NS: Not significant. 
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Significantly detected positively selected lineages are 

highlighted in bold. 

(PDF) 

Table S5. (S5.1-S5.3) Sequences analysed. 

(PDF) 

Figure S1. Maximum-likelihood molecular phylogeny of 
nerve growth factors (NGF). Branches with bootstrap support 
of less than 850 (out of 1000 bootstrap replicates) are 
highlighted in grey. [NFF: "non-front-fanged" advanced snakes; 
Atr: Atractaspidinae ; Sci: Scinciformata]. 
(PDF) 

Figure S2. Maximum-likelihood molecular phylogeny of 
brain-derived neurotrophic factors (BDNF). Branches with 
bootstrap support of less than 850 (out of 1000 bootstrap 
replicates) are highlighted in grey. 
(PDF) 

Figure S3. Maximum-likelihood molecular phylogeny of 
neurotrophin-3 (NT3). Branches with bootstrap support of less 
than 850 (out of 1000 bootstrap replicates) are highlighted in 
grey. 
(PDF) 

Figure S4. Maximum-likelihood phylogeny of 
neurotrophins. Branches with bootstrap support of less than 
850 (out of 1000 bootstrap replicates) are highlighted in grey. 
[NFF: 'non-front-fanged' advanced snakes; Pyt: Pythonidae; 
Mam: Mammals]. 
(PDF) 

Figure S5. Alignment of vertebrate nerve growth factors. 

(PDF) 

Figure S6. Evolutionary fingerprint of nerve growth 
factors (NGF). Estimates of the distribution of synonymous (a) 
and non-synonymous (b) substitution ratesinferred for various 
reptilian and mammalian nerve growth factor (NGF) lineages 
are shown here. The ellipses reflect a Gaussian-approximated 
variance in each individual rate estimate, and coloured pixels 
show the density of the posterior sample of the distribution for a 
given rate. The diagonal line represents the idealized neutral 
evolution regime (w = 1), points above and below the line 
correspond to positive selection (uj>1) and negative selection 
(oj<1 ), respectively. Site model 8 omega (w) along with the 
total number of positively selected sites detected by its Bayes 
Empirical Bayes (BEB) approach are also indicated below. 
(PDF) 

Figure S7. Evolutionary fingerprint of brain-derived 
neurotrophic factors (BDNF). Estimates of the distribution of 



synonymous (a) and non-synonymous (b) substitution rates 
inferred for various reptilian and mammalian brain-derived 
neurotrophic factor (BDNF) lineages are shown here. The 
ellipses reflect a Gaussian-approximated variance in each 
individual rate estimate, and coloured pixels show the density 
of the posterior sample of the distribution for a given rate. The 
diagonal line represents the idealized neutral evolution regime 
(to = 1), points above and below the line correspond to positive 
selection (co>1) and negative selection (w<1), respectively. Site 
model 8 omega (w) along with the total number of positively 
selected sites detected by its Bayes Empirical Bayes (BEB) 
approach are also indicated below. 
(PDF) 

Figure S8. Evolutionary fingerprint of neurotrophin-3 
(NT3). Estimates of the distribution of synonymous (a) and 
non-synonymous (b) substitution rates inferred for various 
reptilian and mammalian neurotrophin-3 (NT-3) lineages are 
shown here. The ellipses reflect a Gaussian-approximated 
variance in each individual rate estimate, and coloured pixels 
show the density of the posterior sample of the distribution for a 
given rate. The diagonal line represents the idealized neutral 
evolution regime (w = 1), points above and below the line 
correspond to positive selection (w>1) and negative selection 
(io<1), respectively. Site model 8 omega (w) along with the 
total number of positively selected sites detected by its Bayes 
Empirical Bayes (BEB) approach is also indicated below. 
(PDF) 

Figure S9. Branch-site REL: Caenophidian nerve growth 
factors (NGF). The hue of each colour indicates strength of 
selection, with primary red corresponding to w > 5, primary 
blue to to = 0 and grey to u>=1. The width of each colour 
component represents the proportion of sites in the 
corresponding class. Thicker branches have been classified as 
undergoing episodic diversifying selection (indicated by arrows) 
by the sequential likelihood ratio test at corrected p < 0.05. 
(PDF) 

File S1. Unpublished sequences. 

(DOCX) 
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